See code here.
#|label: simulation-specifics
## simulation design specifics
# specify the sample sizes
N <- c(50, 150, 500, 1000, 2000, 3000, 4000, 5000, 7500, 10000)
# specify replication number
n <- 500
# specify alpha level
alpha <- 0.01Discovering Cyclic Causal Models
Supplementary Material
1 Introduction
Statistical network models have gained popularity in analyzing multivariate psychological data. These models often interpret network parameters as reflecting causal relationships, making it a form of causal discovery. However, recent research has shown that network models may not perform well as causal discovery tools for discovering acyclic causal structures (DAGs), and alternative methods are preferred for this task (Ryan, Bringmann, and Schuurman 2022).
But, acyclic causal models may not be suitable for some psychological phenomena, such as psychopathologies, which we expect to have cycles or feedback loop relationships between symptoms. While cyclic causal discovery methods have been developed in the computer science literature, they are not widely applied or well understood in empirical practice.
To address this gap, the main paper provides an accessible introduction to the basics of cyclic causal discovery for empirical researchers (Park 2023). It examines three different cyclic causal discovery methods and investigates their performance in typical psychological research contexts through a simulation study.
This supplementary material provides a more detailed analysis of the main simulation study. We omitted these results from the paper due to space constraints. In this supplementary material, we show the estimated PAGs and their associated measures of accuracy and uncertainty for each simulation condition. Additionally, we provide extra results on the running time of the algorithms used in the study.
2 Simulation Conditions
Below, we show the directed cyclic graph (DCG) and the corresponding true ancestral graph for each condition. There are in total eight conditions: model size (\(p = 5\), \(p = 10\)) \(\times\) density (sparse, dense) \(\times\) presence of latent confounder (presence, absence).
2.1 5p_sparse
See code here.
#|label: 5psparse
## ====================
## 5p - sparse
## ====================
# specify B matrix
B5sparse = matrix(c(0, 0, 0, 0, 0,
1, 0, 0.8, 0, 0,
0, 0, 0, 0.9, 0,
0, 0.7, 0, 0, 1.5,
0, 0, 0, 0, 0), 5, 5, byrow = T)
colnames(B5sparse) <- c("X1", "X2", "X3", "X4", "X5")
# specify layout
layout5 = matrix(c(0,1,
0,0,
1,-1,
2,0,
2,1),5,2,byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true5psparse <- qgraph(t(B5sparse), layout=layout5, labels = colnames(B5sparse), theme="colorblind")
title("True cyclic graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)
## Data generating
# equilibrium check
equilibrium_check(B5sparse)
# generate data
simdata_5psparse <- N %>% future_map(function(z) {
replicate(n = n,
expr = gen_dat(B5sparse, N = z),
simplify = FALSE)
}, .options = furrr_options(seed=123))
## True Ancestral Graph
dcg_5psparse <- matrix(c(0,1,0,0,0,
0,0,0,1,0,
0,1,0,0,0,
0,0,1,0,0,
0,0,0,1,0), 5,5,byrow=T)
trueag_5psparse <- true_ancestral(dcg_5psparse, gen_dat(B5sparse), gaussCItest)
dimnames(trueag_5psparse) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(trueag_5psparse)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)
See code here.
# Run CCD algorithm
# ccd_5psparse <- simdata_5psparse %>%
# map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
# )
# mat_5psparse <- ccd_5psparse %>%
# map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_5psparse.RData")
# can we have a more elegant way to combine "graphNEL" plots? ? ? ? ?
# par(mfrow=c(2,5))
# pag_ccd5psparse <- map2(ccd_5psparse, mat_5psparse,
# ~map2(.x, .y, plotPAG)
# )
## Run FCI algorithm
# fci_5psparse <- simdata_5psparse %>%
# map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
# alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # extract amat
# )
load("../data/fci_5psparse.RData")
## Run CCI algorithm
# cci_5psparse <- simdata_5psparse %>%
# map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
# )
load("../data/cci_5psparse.RData")
# par(mfrow=c(2,5))
# pag_cci5psparse <- cci_5psparse %>%
# map_depth(2, ~plotAG(.x))
## evaluation
# CCD
res_ccd5psparse <- mat_5psparse %>%
map_depth(2, ~precision_recall(trueag_5psparse, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_ccd5psparse <- mat_5psparse %>%
map_depth(2, ~uncertainty(.x)) %>% do.call("cbind", .) %>% apply(., 2, unlist) %>% as.data.frame %>%
rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd5psparse <- mat_5psparse %>%
map_depth(2, ~SHD(trueag_5psparse, .x)) %>% do.call("cbind", .) %>% apply(., 2, unlist) %>% as.data.frame %>% rename_with(~ paste0("N = ", N))
# FCI
res_fci5psparse <- fci_5psparse %>%
map_depth(2, ~precision_recall(trueag_5psparse, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_fci5psparse <- fci_5psparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci5psparse <- fci_5psparse %>%
map_depth(2, ~SHD(trueag_5psparse, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# CCI
res_cci5psparse <- cci_5psparse %>%
map_depth(2, ~precision_recall(trueag_5psparse, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_cci5psparse <- cci_5psparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci5psparse <- cci_5psparse %>%
map_depth(2, ~SHD(trueag_5psparse, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))2.2 5p_dense
See code here.
#|label: 5pdense
## ====================
## 5p - dense
## ====================
# specify B matrix
## cyclic case
B5dense = matrix(c(0, 0, 0, 0, 0,
1.4, 0, 0.8, 0, 0,
0, 0, 0, 0.9, 0,
0, 0.7, 0, 0, 1,
1, 0, 0, 0, 0), 5, 5, byrow = T)
colnames(B5dense) <- c("X1", "X2", "X3", "X4", "X5")
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true5pdense <- qgraph(t(B5dense), layout=layout5, labels = colnames(B5dense), theme="colorblind")
title("True cyclic graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)
## Data generating
# equilibrium check
equilibrium_check(B5dense)
# generate data (sample size as specified above)
simdata_5pdense <- N %>% future_map(function(z) {
replicate(n = n,
expr = gen_dat(B5dense, N = z),
simplify = FALSE)
}, .options = furrr_options(seed=123))
# True Ancestral Graph
dcg_5pdense <- matrix(c(0,1,0,0,1,
0,0,0,1,0,
0,1,0,0,0,
0,0,1,0,0,
0,0,0,1,0), 5,5,byrow=T)
trueag_5pdense <- true_ancestral(dcg_5pdense, gen_dat(B5dense), gaussCItest)
dimnames(trueag_5pdense) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(trueag_5pdense)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)
See code here.
## Run CCD algorithm
# ccd_5pdense <- simdata_5pdense %>%
# map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
# )
# mat_5pdense <- ccd_5pdense %>%
# map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
# par(mfrow=c(2,5))
# pag_ccd5pdense <- map2(ccd_5pdense, mat_5pdense,
# ~map2(.x, .y, plotPAG)
# )
load("../data/mat_5pdense.RData")
## Run FCI algorithm
# fci_5pdense <- simdata_5pdense %>%
# map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
# alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # extract amat
# )
load("../data/fci_5pdense.RData")
## Run CCI algorithm
# cci_5pdense <- simdata_5pdense %>%
# map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
# )
load("../data/cci_5pdense.RData")
## evaluation
# CCD
res_ccd5pdense <- mat_5pdense %>%
map_depth(2, ~precision_recall(trueag_5pdense, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_ccd5pdense <- mat_5pdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd5pdense <- mat_5pdense %>%
map_depth(2, ~SHD(trueag_5pdense, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# FCI
res_fci5pdense <- fci_5pdense %>%
map_depth(2, ~precision_recall(trueag_5pdense, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_fci5pdense <- fci_5pdense%>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci5pdense <- fci_5pdense %>%
map_depth(2, ~SHD(trueag_5pdense, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# CCI
res_cci5pdense <- cci_5pdense %>%
map_depth(2, ~precision_recall(trueag_5pdense, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_cci5pdense <- cci_5pdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci5pdense <- cci_5pdense %>%
map_depth(2, ~SHD(trueag_5pdense, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))2.3 10p_sparse
See code here.
## ====================
## 10p - sparse
## ====================
B10sparse = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0.8, 0, 0, 0, 0, 0, 0, 0,
0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0.7, 0, 0, 0.9, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0.8, 0, 0.5, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.8, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 10, 10, byrow = T)
dimnames(B10sparse) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
# specify layout
layout10 = matrix(c(0,1,
2,1,
1,0,
2,-1,
3,0,
4, -1,
5, 0,
6, -1,
4, 1,
7, 1),10,2,byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true10psparse <- qgraph(t(B10sparse), layout = layout10, labels = colnames(B10sparse), theme="colorblind")
title("True cyclic graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)
## Data generating
# equilibrium check
equilibrium_check(B10sparse)
# generate data (sample size as specified above)
simdata_10psparse <- N %>% future_map(function(z) {
replicate(n = n,
expr = gen_dat(B10sparse, N = z),
simplify = FALSE)
}, .options = furrr_options(seed=123))
## True Ancestral Graph
dcg_10psparse <- matrix(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0), 10, 10, byrow = T)
trueag_10psparse <- true_ancestral(dcg_10psparse, gen_dat(B10sparse), gaussCItest)
dimnames(trueag_10psparse) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(trueag_10psparse)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)
See code here.
## Run CCD algorithm
# ccd_10psparse <- simdata_10psparse %>%
# map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
# )
#
# mat_10psparse <- ccd_10psparse %>%
# map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_10psparse.RData")
## Run FCI algorithm
# fci_10psparse <- simdata_10psparse %>%
# map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
# alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
# )
load("../data/fci_10psparse.RData")
## Run CCI algorithm
# cci_10psparse <- simdata_10psparse %>%
# map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
# )
load("../data/cci_10psparse.RData")
## evaluation
# CCD
res_ccd10psparse <- mat_10psparse %>%
map_depth(2,
~precision_recall(trueag_10psparse, .x)) %>%
do.call("cbind", .) %>% t() %>%
apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_ccd10psparse <- mat_10psparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd10psparse <- mat_10psparse %>%
map_depth(2, ~SHD(trueag_10psparse, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# FCI
res_fci10psparse <- fci_10psparse %>%
map_depth(2, ~precision_recall(trueag_10psparse, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_fci10psparse <- fci_10psparse%>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci10psparse <- fci_10psparse %>%
map_depth(2, ~SHD(trueag_10psparse, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# CCI
res_cci10psparse <- cci_10psparse %>%
map_depth(2, ~precision_recall(trueag_10psparse, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_cci10psparse <- cci_10psparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci10psparse <- cci_10psparse %>%
map_depth(2, ~SHD(trueag_10psparse, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))2.4 10p_dense
See code here.
## ====================
## 10p - dense
## ====================
B10dense = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.3, 0, 0.8, 0, 0, 0, 0, 0, 0, 0,
0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1.1, 0, 0, 0.9, 0, 0, 0, 0,
0, 0.4, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0.9, 0, 0.5, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.8, 1,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0.4,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.7,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 10, 10, byrow = T)
colnames(B10dense) <- c(paste("X", 1:10, sep=""))
# specify layout
layout10 = matrix(c(0,1,
2,1,
1,0,
2,-1,
3,0,
4, -1,
5, 0,
6, -1,
4, 1,
7, 1),10,2,byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true10pdense <- qgraph(t(B10dense), layout = layout10, labels = colnames(B10dense), theme="colorblind")
title("True cyclic graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)
## Data generating
# equilibrium check
equilibrium_check(B10dense)
# generate data (sample size as specified above)
simdata_10pdense <- N %>% future_map(function(z) {
replicate(n = n,
expr = gen_dat(B10dense, N = z),
simplify = FALSE)
}, .options = furrr_options(seed=123))
## True Ancestral Graph
dcg_10pdense <- matrix(c(0, 1, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 1, 1, 0), 10, 10, byrow = T)
trueag_10pdense <- true_ancestral(dcg_10pdense, gen_dat(B10dense), gaussCItest)
dimnames(trueag_10pdense) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(trueag_10pdense)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)
See code here.
## Run CCD algorithm
# ccd_10pdense <- simdata_10pdense %>%
# map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
# )
# mat_10pdense <- ccd_10pdense %>%
# map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_10pdense.RData")
## Run FCI algorithm
# fci_10pdense <- simdata_10pdense %>%
# map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
# alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
# )
load("../data/fci_10pdense.RData")
## Run CCI algorithm
# cci_10pdense <- simdata_10pdense %>%
# map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
# )
load("../data/cci_10pdense.RData")
## evaluation
# CCD
res_ccd10pdense <- mat_10pdense %>%
map_depth(2, ~precision_recall(trueag_10pdense, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_ccd10pdense <- mat_10pdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd10pdense <- mat_10pdense %>%
map_depth(2, ~SHD(trueag_10pdense, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# FCI
res_fci10pdense <- fci_10pdense %>%
map_depth(2, ~precision_recall(trueag_10pdense, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_fci10pdense <- fci_10pdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci10pdense <- fci_10pdense %>%
map_depth(2, ~SHD(trueag_10pdense, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# CCI
res_cci10pdense <- cci_10pdense %>%
map_depth(2, ~precision_recall(trueag_10pdense, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_cci10pdense <- cci_10pdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci10pdense <- cci_10pdense %>%
map_depth(2, ~SHD(trueag_10pdense, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))2.5 5p_LVsparse
See code here.
## ====================
## 5p with LV sparse
## ====================
# specify B matrix
B5_lvsparse = matrix(c(0, 0, 0, 0, 0, 1,
0, 0, 0.4, 0, 0, 1,
0, 0, 0, 0.5, 0,0,
0, 0.7, 0, 0, 1.5,0,
0, 0, 0, 0, 0,0,
0,0,0,0,0,0), 6, 6, byrow = T)
colnames(B5_lvsparse) <- c("X1", "X2", "X3", "X4", "X5", "L1")
# specify layout
layout5_lv = matrix(c(0,1,
0,0,
1,-1,
2,0,
2,1,
-1, 0.5),6,2,byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true5p_lvsparse <- qgraph(t(B5_lvsparse), layout=layout5_lv, labels = colnames(B5_lvsparse), theme="colorblind")
title("True cyclic graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)
## Data generating
# equilibrium check
equilibrium_check(B5_lvsparse)
# generate data (sample size as specified above)
simdata_5pLVsparse <- N %>% future_map(function(z) {
replicate(n = n,
expr = gen_dat(B5_lvsparse, N = z)[,-6],
simplify = FALSE)
}, .options = furrr_options(seed=123))
## True Ancestral Graph
# [i,j] = [j,i] = 2: a LV exists between i and j
dcg_5psparseLV <- matrix(c(0, 2, 0, 0, 0,
2, 0, 0, 1, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0), 5, 5, byrow = T)
trueag_5psparseLV <- true_ancestral(dcg_5psparseLV, gen_dat(B5_lvsparse), gaussCItest)
dimnames(trueag_5psparseLV) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(trueag_5psparseLV)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)
See code here.
## Run CCD algorithm
# ccd_5pLVsparse <- simdata_5pLVsparse %>%
# map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
# )
# mat_5pLVsparse <- ccd_5pLVsparse %>%
# map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_5pLVsparse.RData")
## Run FCI algorithm
# fci_5pLVsparse <- simdata_5pLVsparse %>%
# map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
# alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
# )
load("../data/fci_5pLVsparse.RData")
## Run CCI algorithm
# cci_5pLVsparse <- simdata_5pLVsparse %>%
# map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
# )
load("../data/cci_5pLVsparse.RData")
## evaluation
# CCD
res_ccd5pLVsparse <- mat_5pLVsparse %>%
map_depth(2, ~precision_recall(trueag_5psparseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_ccd5pLVsparse <- mat_5pLVsparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd5pLVsparse <- mat_5pLVsparse %>%
map_depth(2, ~SHD(trueag_5psparseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# FCI
res_fci5pLVsparse <- fci_5pLVsparse %>%
map_depth(2, ~precision_recall(trueag_5psparseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_fci5pLVsparse <- fci_5pLVsparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci5pLVsparse <- fci_5pLVsparse %>%
map_depth(2, ~SHD(trueag_5psparseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# CCI
res_cci5pLVsparse <- cci_5pLVsparse %>%
map_depth(2, ~precision_recall(trueag_5psparseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_cci5pLVsparse <- cci_5pLVsparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci5pLVsparse <- cci_5pLVsparse %>%
map_depth(2, ~SHD(trueag_5psparseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))2.6 5p_LVdense
See code here.
## ====================
## 5p with LV dense
## ====================
# specify B matrix
B5_lvdense = matrix(c(0, 0, 0, 0, 0, 1,
0, 0, 0.4, 0, 0, 1,
0, 0, 0, 0.5, 0,0,
0, 0.7, 0, 0, 0.5, 0,
0.6, 0, 0, 0, 0,0,
0,0,0,0,0,0), 6, 6, byrow = T)
colnames(B5_lvdense) <- c("X1", "X2", "X3", "X4", "X5", "L1")
# specify layout
layout5_lv = matrix(c(0,1,
0,0,
1,-1,
2,0,
2,1,
-1, 0.5),6,2,byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true5p_lvdense <- qgraph(t(B5_lvdense), layout=layout5_lv, labels = colnames(B5_lvdense), theme="colorblind")
title("True cyclic graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)
## Data generating
# equilibrium check
equilibrium_check(B5_lvdense)
# generate data (sample size as specified above)
simdata_5pLVdense <- N %>% future_map(function(z) {
replicate(n = n,
expr = gen_dat(B5_lvdense, N = z)[,-6],
simplify = FALSE)
}, .options = furrr_options(seed=123))
## True Ancestral Graph
# [i,j] = [j,i] = 2: a LV exists between i and j
dcg_5pdenseLV <- matrix(c(0, 2, 0, 0, 1,
2, 0, 0, 1, 0,
0, 1, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 1, 0), 5, 5, byrow = T)
trueag_5pdenseLV <- true_ancestral(dcg_5pdenseLV, gen_dat(B5_lvdense), gaussCItest)
dimnames(trueag_5pdenseLV) <- list(paste("X", 1:5, sep=""), paste("X", 1:5, sep=""))
plotAG(trueag_5pdenseLV)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)
See code here.
## Run CCD algorithm
# ccd_5pLVdense <- simdata_5pLVdense %>%
# map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
# )
# mat_5pLVdense <- ccd_5pLVdense %>%
# map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_5pLVdense.RData")
## Run FCI algorithm
# fci_5pLVdense <- simdata_5pLVdense %>%
# map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
# alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
# )
#
load("../data/fci_5pLVdense.RData")
## Run CCI algorithm
# cci_5pLVdense <- simdata_5pLVdense %>%
# map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
# )
#
load("../data/cci_5pLVdense.RData")
## evaluation
# CCD
res_ccd5pLVdense <- mat_5pLVdense %>%
map_depth(2, ~precision_recall(trueag_5pdenseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_ccd5pLVdense <- mat_5pLVdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd5pLVdense <- mat_5pLVdense %>%
map_depth(2, ~SHD(trueag_5pdenseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# FCI
res_fci5pLVdense <- fci_5pLVdense %>%
map_depth(2, ~precision_recall(trueag_5pdenseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_fci5pLVdense <- fci_5pLVdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci5pLVdense <- fci_5pLVdense %>%
map_depth(2, ~SHD(trueag_5pdenseLV , .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# CCI
res_cci5pLVdense <- cci_5pLVdense %>%
map_depth(2, ~precision_recall(trueag_5pdenseLV , .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_cci5pLVdense <- cci_5pLVdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci5pLVdense <- cci_5pLVdense %>%
map_depth(2, ~SHD(trueag_5pdenseLV , .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))2.7 10p_LVsparse
See code here.
## ====================
## 10p with LV sparse
## ====================
# specify B matrix
## with 2 LVs
B10_lvsparse = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0.8, 0, 0, 0, 0, 0, 0, 0, 0, 0.7,
0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0.7, 0, 0, 0.9, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.7,
0, 0, 0, 0, 0.8, 0, 0.5, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0.8, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.8, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 12, 12, byrow = T)
colnames(B10_lvsparse) <- c(paste("X", 1:10, sep=""), "L1", "L2")
# specify layout
layout10LV2 = matrix(c(0, 1,
2, 1,
1, 0,
2, -1,
3, 0,
4, -1,
5, 0,
6, -1,
4, 1,
7, 1,
8, 0,
3, 2), 12, 2, byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true10pLVsparse <- qgraph(t(B10_lvsparse), layout = layout10LV2, labels = colnames(B10_lvsparse), theme="colorblind")
title("True cyclic graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)
## Data generating
# equilibrium check
equilibrium_check(B10_lvsparse)
# generate data (sample size as specified above)
simdata_10pLVsparse <- N %>% future_map(function(z) {
replicate(n = n,
expr = gen_dat(B10_lvsparse, N = z)[,-c(11,12)],
simplify = FALSE)
}, .options = furrr_options(seed=123))
## True Ancestral Graph
# [i,j] = [j,i] = 2: a LV exists between i and j
dcg_10psparseLV <- matrix(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 2, 0, 0, 0, 0, 0,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 2, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 2,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 2, 1, 0), 10, 10, byrow = T)
trueag_10psparseLV <- true_ancestral(dcg_10psparseLV, gen_dat(B10_lvsparse), gaussCItest)
dimnames(trueag_10psparseLV) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(trueag_10psparseLV)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)
See code here.
## Run CCD algorithm
# ccd_10pLVsparse <- simdata_10pLVsparse %>%
# map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
# )
# mat_10pLVsparse <- ccd_10pLVsparse %>%
# map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_10pLVsparse.RData")
## Run FCI algorithm
# fci_10pLVsparse <- simdata_10pLVsparse %>%
# map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
# alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
# )
load("../data/fci_10pLVsparse.RData")
## Run CCI algorithm
# cci_10pLVsparse <- simdata_10pLVsparse %>%
# map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
# )
load("../data/cci_10pLVsparse.RData")
## evaluation
# CCD
res_ccd10pLVsparse <- mat_10pLVsparse %>%
map_depth(2, ~precision_recall(trueag_10psparseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_ccd10pLVsparse <- mat_10pLVsparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd10pLVsparse <- mat_10pLVsparse %>%
map_depth(2, ~SHD(trueag_10psparseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# FCI
res_fci10pLVsparse <- fci_10pLVsparse %>%
map_depth(2, ~precision_recall(trueag_10psparseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_fci10pLVsparse <- fci_10pLVsparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci10pLVsparse <- fci_10pLVsparse %>%
map_depth(2, ~SHD(trueag_10psparseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# CCI
res_cci10pLVsparse <- cci_10pLVsparse %>%
map_depth(2, ~precision_recall(trueag_10psparseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_cci10pLVsparse <- cci_10pLVsparse %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci10pLVsparse <- cci_10pLVsparse %>%
map_depth(2, ~SHD(trueag_10psparseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))2.8 10p_LVdense
See code here.
## ====================
## 10p with LV dense
## ====================
# specify B matrix
B10_lvdense = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0.5, 0, 1.4, 0, 0, 0, 0, 0, 0, 0, 0, 1.1,
0.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0.9, 0, 0, 0, 0, 0, 0,
0, 0, 1.2, 0.8, 0, 0, 0, 0, 0, 0, 0, 0.7,
0, 0, 0, 0, 0.8, 0, 0.5, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1.4, 0, 0.4, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.6, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.6, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 12, 12, byrow = T)
colnames(B10_lvdense) <- c(paste("X", 1:10, sep=""), "L1", "L2")
# specify layout
layout10LV2 = matrix(c(0, 1,
2, 1,
1, 0,
2, -1,
3, 0,
4, -1,
5, 0,
6, -1,
4, 1,
7, 1,
8, 0,
3, 2), 12, 2, byrow = T)
par(oma=c(0, 0, 4, 0), mfrow=c(1,2))
true10pLVdense <- qgraph(t(B10_lvdense), layout = layout10LV2, labels = colnames(B10_lvdense), theme="colorblind")
title("True cyclic graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.18)
## Data generating
# equilibrium check
equilibrium_check(B10_lvdense)
# generate data (sample size as specified above)
simdata_10pLVdense <- N %>% future_map(function(z) {
replicate(n = n,
expr = gen_dat(B10_lvdense, N = z)[,-c(11,12)],
simplify = FALSE)
}, .options = furrr_options(seed=123))
## True Ancestral Graph
# [i,j] = [j,i] = 2: a LV exists between i and j
dcg_10pdenseLV <- matrix(c(0, 1, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 2, 0, 0, 0, 0, 0,
0, 1, 0, 1, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 2, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 1, 0, 2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 2,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 2, 2, 1, 0), 10, 10, byrow = T)
trueag_10pdenseLV <- true_ancestral(dcg_10pdenseLV, gen_dat(B10_lvdense), gaussCItest)
dimnames(trueag_10pdenseLV) <- list(paste("X", 1:10, sep=""), paste("X", 1:10, sep=""))
plotAG(trueag_10pdenseLV)
title(main = "True ancestral graph", font.main = 1, cex.main = 1.2, line = 2, outer=TRUE, adj = 0.86)
See code here.
## Run CCD algorithm
# ccd_10pLVdense <- simdata_10pLVdense %>%
# map_depth(2, ~ ccdKP(df = .x, dataType = "continuous", alpha = alpha)
# )
# mat_10pLVdense <- ccd_10pLVdense %>%
# map_depth(2, ~CreateAdjMat(.x, length(.x$nodes)))
load("../data/mat_10pLVdense.RData")
## Run FCI algorithm
# fci_10pLVdense <- simdata_10pLVdense %>%
# map_depth(2, ~fci(list(C = cor(.x), n = nrow(.x)), indepTest=gaussCItest,
# alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(.x)) %>% .@amat # exxtract amat
# )
load("../data/fci_10pLVdense.RData")
## Run CCI algorithm
# cci_10pLVdense <- simdata_10pLVdense %>%
# map_depth(2, ~cci(list(C = cor(.x), n = nrow(.x)), gaussCItest, alpha=alpha, labels = colnames(.x), p = ncol(.x)) %>% .$maag # convert some logical matrix (0, 1 only) to a numeric matrix while keeping a matrix format (lost the row names but they are not needed)
# )
load("../data/cci_10pLVdense.RData")
## evaluation
# CCD
res_ccd10pLVdense <- mat_10pLVdense %>%
map_depth(2, ~precision_recall(trueag_10pdenseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_ccd10pLVdense <- mat_10pLVdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_ccd10pLVdense <- mat_10pLVdense %>%
map_depth(2, ~SHD(trueag_10pdenseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# FCI
res_fci10pLVdense <- fci_10pLVdense %>%
map_depth(2, ~precision_recall(trueag_10pdenseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_fci10pLVdense <- fci_10pLVdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_fci10pLVdense <- fci_10pLVdense %>%
map_depth(2, ~SHD(trueag_10pdenseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# CCI
res_cci10pLVdense <- cci_10pLVdense %>%
map_depth(2, ~precision_recall(trueag_10pdenseLV, .x)) %>%
do.call("cbind", .) %>% t() %>% apply(., 2, unlist) %>% as.data.frame()
# UNCERTAINTY
uncer_cci10pLVdense <- cci_10pLVdense %>%
map_depth(2, ~uncertainty(.x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))
# SHD
SHD_cci10pLVdense <- cci_10pLVdense %>%
map_depth(2, ~SHD(trueag_10pdenseLV, .x)) %>%
do.call("cbind", .) %>% apply(., 2, unlist) %>%
as.data.frame %>% rename_with(~ paste0("N = ", N))3 Performance Evaluation per Condition
Below, we present the most frequently estimated PAG for each sample size (N) from each algorithm. The graphs are obtained by picking the most frequent type of edge-endpoint produced by algorithms out of 500 simulations. Additionally, we present the matrix plots of the correct estimation and uncertainty for each condition. The darker the color (blue/red) is, the higher the rate of correct estimation (blue) or uncertainty (red).
See code here.
par(mfrow=c(2,5))
vec <- list("CCD-5p-sparse" = mat_5psparse, "FCI-5p-sparse" = fci_5psparse, "CCI-5p-sparse" = cci_5psparse)
graphs_5psparse <- list()
for(i in seq_along(vec)){
## high-freq graphs
graphs_5psparse[[i]] <- vec[[i]] %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
## correct prop plots
vec[[i]] %>%
imap(
~ prop_correct(.x, trueag_5psparse, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob(glue::glue("Correct Proportion {names(vec[[i]])}"), face = "bold", size = 18, family = "Palatino"))
## uncertainty prop plots
vec[[i]] %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob(glue::glue("Uncertainty {names(vec)[i]}"), face = "bold", size = 18, family = "Palatino"))
}
corprop_plots_5psparse <- list()
# prop correct
for(i in seq_along(vec)){
corprop_plots_5psparse[[i]] <- vec[[i]] %>%
imap(
~ prop_correct(.x, trueag_5psparse, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob(glue::glue("Correct Proportion {names(vec[[i]])}"), face = "bold", size = 18, family = "Palatino"))
}
ucprop_plots_5psparse <- list()
# prop uncertain
for(i in seq_along(vec)){
ucprop_plots_5psparse[[i]] <- vec[[i]] %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob(glue::glue("Uncertainty {names(vec)[i]}"), face = "bold", size = 18, family = "Palatino"))
}See code here.
## CCD 5p sparse case
# high frequency
par(mfrow=c(2,5))
mat_5psparse %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
mat_5psparse %>%
imap(
~ prop_correct(.x, trueag_5psparse, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-5p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
mat_5psparse %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-5p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## FCI 5p sparse case
# high frequency
par(mfrow=c(2,5))
fci_5psparse %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
fci_5psparse %>%
imap(
~ prop_correct(.x, trueag_5psparse, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
fci_5psparse %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))
cci_5psparse %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
cci_5psparse %>%
imap(
~ prop_correct(.x, trueag_5psparse, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 6, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
cci_5psparse %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCD 5p dense case
# high frequency
par(mfrow=c(2,5))
mat_5pdense %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
mat_5pdense %>%
imap(
~ prop_correct(.x, trueag_5pdense, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 6, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-5p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
mat_5pdense %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-5p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
## FCI 5p dense case
# high frequency
par(mfrow=c(2,5))
fci_5pdense %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
fci_5pdense %>%
imap(
~ prop_correct(.x, trueag_5pdense, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-5p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
fci_5pdense %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-5p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCI 5p dense case
# high frequency
par(mfrow=c(2,5))
cci_5pdense %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
cci_5pdense %>%
imap(
~ prop_correct(.x, trueag_5pdense, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-5p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
cci_5pdense %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-5p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCD 10p sparse case
# high frequency
par(mfrow=c(2,5))
mat_10psparse %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
mat_10psparse %>%
imap(
~ prop_correct(.x, trueag_10psparse, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-10p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
mat_10psparse %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-10p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## FCI 10p sparse case
# high frequency
par(mfrow=c(2,5))
fci_10psparse %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
fci_10psparse %>%
imap(
~ prop_correct(.x, trueag_10psparse, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
fci_10psparse %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCI 10p sparse case
# high frequency
par(mfrow=c(2,5))
cci_10psparse %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
cci_10psparse %>%
imap(
~ prop_correct(.x, trueag_10psparse, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
cci_10psparse %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCD 10p dense case
# high frequency
par(mfrow=c(2,5))
mat_10pdense %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
mat_10pdense %>%
imap(
~ prop_correct(.x, trueag_10pdense, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 6, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-10p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
mat_10pdense %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-10p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
## FCI 10p dense case
# high frequency
par(mfrow=c(2,5))
fci_10pdense %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
fci_10pdense %>%
imap(
~ prop_correct(.x, trueag_10psparse, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-10p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
fci_10pdense %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-10p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCI 10p dense case
# high frequency
par(mfrow=c(2,5))
cci_10pdense %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
cci_10pdense %>%
imap(
~ prop_correct(.x, trueag_10psparse, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-10p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
cci_10pdense %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 6, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-10p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCD 5p sparse LV case
# high frequency
par(mfrow=c(2,5))
mat_5pLVsparse %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
mat_5pLVsparse %>%
imap(
~ prop_correct(.x, trueag_5psparseLV, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-5p-LV-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
mat_5pLVsparse %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-5p-LV-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# high frequency
par(mfrow=c(2,5))
fci_5pLVsparse %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
fci_5pLVsparse %>%
imap(
~ prop_correct(.x, trueag_5psparseLV, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
fci_5pLVsparse %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-5p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))
cci_5pLVsparse %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
cci_5pLVsparse %>%
imap(
~ prop_correct(.x, trueag_5psparseLV, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-5p-LV-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
cci_5pLVsparse %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-5p-LV-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCD 5p dense LV case
# high frequency
par(mfrow=c(2,5))
mat_5pLVdense %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
mat_5pLVdense %>%
imap(
~ prop_correct(.x, trueag_5pdenseLV, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-5p-LV-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
mat_5pLVdense %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-5p-LV-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# high frequency
par(mfrow=c(2,5))
fci_5pLVdense %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
fci_5pLVdense %>%
imap(
~ prop_correct(.x, trueag_5pdenseLV, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-5p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
fci_5pLVdense %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-5p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))
cci_5pLVdense %>%
map(~high_freq(.x, p = 5) %>%
plotAG)
See code here.
# prop correct
cci_5pLVdense %>%
imap(
~ prop_correct(.x, trueag_5pdenseLV, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-5p-LV-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
cci_5pLVdense %>%
imap(
~ prop_uncertain(.x, p = 5) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-5p-LV-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCD 10p sparse LV case
# high frequency
par(mfrow=c(2,5))
mat_10pLVsparse %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
mat_10pLVsparse %>%
imap(
~ prop_correct(.x, trueag_10psparseLV, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-10p-LV-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
mat_10pLVsparse %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-10p-LV-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# high frequency
par(mfrow=c(2,5))
fci_10pLVsparse %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
fci_10pLVsparse %>%
imap(
~ prop_correct(.x, trueag_10psparseLV, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
fci_10pLVsparse %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-10p-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))
cci_10pLVsparse %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
cci_10pLVsparse %>%
imap(
~ prop_correct(.x, trueag_10psparseLV, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-10p-LV-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
cci_10pLVsparse %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-10p-LV-sparse", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCD 10p dense LV case
# high frequency
par(mfrow=c(2,5))
mat_10pLVdense %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
mat_10pLVdense %>%
imap(
~ prop_correct(.x, trueag_10pdenseLV, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCD-10p-LV-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
mat_10pLVdense %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCD-10p-LV-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# high frequency
par(mfrow=c(2,5))
fci_10pLVdense %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
fci_10pLVdense %>%
imap(
~ prop_correct(.x, trueag_10pdenseLV, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion FCI-10p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
fci_10pLVdense %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty FCI-10p-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
## CCI 5p sparse case
# high frequency
par(mfrow=c(2,5))
cci_10pLVdense %>%
map(~high_freq(.x, p = 10) %>%
plotAG)
See code here.
# prop correct
cci_10pLVdense %>%
imap(
~ prop_correct(.x, trueag_10pdenseLV, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="blue") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Correct Proportion CCI-10p-LV-dense", face = "bold", size = 18, family = "Palatino"))
See code here.
# prop uncertain
cci_10pLVdense %>%
imap(
~ prop_uncertain(.x, p = 10) %>%
# long format
reshape2::melt() %>%
ggplot(aes(x = Var2, y = Var1, fill = value, label= value)) +
geom_tile() +
geom_text() +
# reverse factor level
scale_y_discrete(limits=rev) +
scale_fill_gradient(low="grey90", high="red") +
labs(x = "", y="", title = glue::glue("N = {N[.y]}"))
) %>%
ggpubr::ggarrange(plotlist = .,
ncol = 5, nrow = 2, common.legend = TRUE, legend = "bottom") %>%
ggpubr::annotate_figure(top = text_grob("Uncertainty CCI-10p-LV-dense", face = "bold", size = 18, family = "Palatino"))
4 Overall Performance Evalulation
Here we summarize the performance of the algorithms across all conditions using the evaluation metrics: structural Hamming distance (SHD), precision, recall, and uncertainty. Each point represents the average value of each metric and shaded area represents interquartile range (IQR).
See code here.
## ============================
## Create neat dataframe
## ============================
## Compute average precision & recall and corresponding sd for each condition
pre_rec <- list(
# put all the results together in a list
res_ccd5psparse, res_fci5psparse, res_cci5psparse, res_ccd10psparse, res_fci10psparse, res_cci10psparse, res_ccd5pdense, res_fci5pdense, res_cci5pdense, res_ccd10pdense, res_fci10pdense, res_cci10pdense, res_ccd5pLVsparse, res_fci5pLVsparse, res_cci5pLVsparse, res_ccd5pLVdense, res_fci5pLVdense, res_cci5pLVdense, res_ccd10pLVsparse, res_fci10pLVsparse, res_cci10pLVsparse, res_ccd10pdense, res_fci10pLVdense, res_cci10pLVdense
) %>%
# transpose df
map(~ sjmisc::rotate_df(.x) %>%
# add sample size (N) info
rename_with(~paste0(.x, "N = ", rep(N, each=8))) %>%
# think about how to deal with NAs or do I want to define sth. else instead of NAs.
# na.omit(.x) %>%
# get the average and sd
dplyr::summarise(across(everything(.), list(mean = ~mean(., na.rm=T), sd = ~sd(., na.rm=T))))) %>%
bind_rows() %>%
mutate(algorithm = rep(c("CCD", "FCI", "CCI"), 8),
condition = rep(c("5p_sparse", "10p_sparse", "5p_dense", "10p_dense", "5p_LVsparse", "5p_LVdense", "10p_LVsparse", "10p_LVdense"), each=3),
netsize = stringr::str_split(condition, "_", simplify=T)[,1],
latentvar = ifelse(stringr::str_detect(condition, "LV")==TRUE, "with LC", "without LC"),
densities = stringr::str_remove(stringr::str_split(condition, "_", simplify=T)[,2], "LV")
) %>%
# brings the algorithm and condition names first
relocate(where(is.character), .before = where(is.numeric)) %>%
# convert it to a long format
tidyr::pivot_longer(!c(algorithm, condition, netsize, latentvar, densities), names_to = "metric", values_to = "value") %>%
# Add sample size column (N) & clean up the column name
mutate(N = stringr::str_extract(metric, "(?<=[N =])\\d+"),
metric = stringr::str_replace_all(metric, "[0-9.]+|[N =]", ""))
## Compute average uncertainty rate and corresponding sd for each condition
uncertainties <- bind_rows(
# bind all results from each condition
"CCD_5p-sparse" = uncer_ccd5psparse, "FCI_5p-sparse" = uncer_fci5psparse, "CCI_5p-sparse"=uncer_cci5psparse, "CCD_10p-sparse"=uncer_ccd10psparse, "FCI_10p-sparse" = uncer_fci10psparse, "CCI_10p-sparse" = uncer_cci10psparse, "CCD_5p-dense"=uncer_ccd5pdense, "FCI_5p-dense"=uncer_fci5pdense, "CCI_5p-dense"=uncer_cci5pdense, "CCD_10p-dense"=uncer_ccd10pdense, "FCI_10p-dense"=uncer_fci10pdense, "CCI_10p-dense"=uncer_cci10pdense, "CCD_5p-LVsparse"=uncer_ccd5pLVsparse, "FCI_5p-LVsparse"=uncer_fci5pLVsparse, "CCI_5p-LVsparse"=uncer_cci5pLVsparse, "CCD_10p-LVsparse"=uncer_ccd10pLVsparse, "FCI_10p-LVsparse"=uncer_fci10pLVsparse, "CCI_10p-LVsparse"=uncer_cci10pLVsparse,
"CCD_5p-LVdense"=uncer_ccd5pLVdense, "FCI_5p-LVdense"=uncer_fci5pLVdense, "CCI_5p-LVdense"=uncer_cci5pLVdense, "CCD_10p-LVdense"=uncer_ccd10pLVdense, "FCI_10p-LVdense"=uncer_fci10pLVdense, "CCI_10p-LVdense"=uncer_cci10pLVdense, .id="id"
) %>%
group_by(id) %>%
# get the average and sd
summarise_all(list(means = mean, sds = sd)) %>%
mutate(algorithm = stringr::str_split(id, "_", simplify = T)[,1],
condition = stringr::str_split(id, "_", simplify = T)[,2],
netsize = stringr::str_split(condition, "-", simplify=T)[,1],
latentvar = ifelse(stringr::str_detect(condition, "LV")==TRUE, "with LC", "without LC"),
densities = stringr::str_remove(stringr::str_split(condition, "-", simplify=T)[,2], "LV")
) %>%
tidyr::pivot_longer(!c(algorithm, condition, id, netsize, latentvar, densities), names_to = "name", values_to = "value") %>%
mutate(N = stringr::str_extract(stringr::str_split(name, "_", simplify = T)[,1], "(\\d)+"),
statistics = stringr::str_split(name, "_", simplify = T)[,2]) %>%
dplyr::select(-id, -name) %>% relocate(where(is.character), .before = where(is.numeric))
## Compute average SHD values and corresponding sd for each condition
SHDs <- bind_rows(
# bind all results from each condition
"CCD_5p-sparse" = SHD_ccd5psparse, "FCI_5p-sparse" = SHD_fci5psparse, "CCI_5p-sparse"=SHD_cci5psparse, "CCD_10p-sparse"= SHD_ccd10psparse, "FCI_10p-sparse" = SHD_fci10psparse, "CCI_10p-sparse" = SHD_cci10psparse, "CCD_5p-dense"= SHD_ccd5pdense, "FCI_5p-dense"=SHD_fci5pdense, "CCI_5p-dense"=SHD_cci5pdense, "CCD_10p-dense"= SHD_ccd10pdense, "FCI_10p-dense"=SHD_fci10pdense, "CCI_10p-dense"=SHD_cci10pdense, "CCD_5p-LVsparse"=SHD_ccd5pLVsparse, "FCI_5p-LVsparse"=SHD_fci5pLVsparse, "CCI_5p-LVsparse"=SHD_cci5pLVsparse, "CCD_10p-LVsparse"=SHD_ccd10pLVsparse, "FCI_10p-LVsparse"=SHD_fci10pLVsparse, "CCI_10p-LVsparse"=SHD_cci10pLVsparse,
"CCD_5p-LVdense"=SHD_ccd5pLVdense, "FCI_5p-LVdense"=SHD_fci5pLVdense, "CCI_5p-LVdense"=SHD_cci5pLVdense, "CCD_10p-LVdense"=SHD_ccd10pLVdense, "FCI_10p-LVdense"=SHD_fci10pLVdense, "CCI_10p-LVdense"=SHD_cci10pLVdense, .id="id"
) %>%
group_by(id) %>%
# get the average and sd
summarise_all(list(means = mean, sds = sd)) %>%
mutate(algorithm = stringr::str_split(id, "_", simplify = T)[,1],
condition = stringr::str_split(id, "_", simplify = T)[,2],
netsize = stringr::str_split(condition, "-", simplify=T)[,1],
latentvar = ifelse(stringr::str_detect(condition, "LV")==TRUE, "with LC", "without LC"),
densities = stringr::str_remove(stringr::str_split(condition, "-", simplify=T)[,2], "LV")
) %>%
tidyr::pivot_longer(!c(algorithm, condition, id, netsize, latentvar, densities), names_to = "name", values_to = "value") %>%
mutate(N = stringr::str_extract(stringr::str_split(name, "_", simplify = T)[,1], "(\\d)+"),
statistics = stringr::str_split(name, "_", simplify = T)[,2]) %>%
dplyr::select(-id, -name) %>% relocate(where(is.character), .before = where(is.numeric))
## ============================
## Create figures
## ============================
## specify the common figure theme
MyTheme <- theme(plot.title = element_text(face = "bold", family = "Palatino", size = 15, hjust=0.5),
plot.subtitle = element_text(face = "italic", family = "Palatino", size = 15, hjust=0.5),
axis.text=element_text(face = "bold",family = "Palatino", size = 11),
axis.text.x = element_text(angle = 45, hjust = 1.2, vjust =1.2),
axis.title = element_text(face = "bold",family = "Palatino", size = 12),
legend.text = element_text(face = "bold", family = "Palatino", size = 12),
legend.position="bottom",
strip.text = element_text(face="bold", size=13, family = "Palatino"),
strip.background = element_rect(fill="#f0f0f0", linetype = "solid", color="gray"),
strip.placement = "outside",
panel.border = element_rect(color = "#DCDCDC", fill = NA),
panel.spacing.y = unit(4, "mm")
)
## SHD figure
SHDs %>%
tidyr::pivot_wider(names_from = statistics, values_from=value) %>%
ggplot(aes(x= as.numeric(N), y=means, group = algorithm, colour = algorithm, fill = algorithm)) +
# add line plots
geom_line(aes(group = algorithm)) +
# add scattered points
geom_point(size=1) +
# add interquartile range (IQR)
geom_ribbon(aes(ymin=means+qnorm(0.25)*sds, ymax=means+qnorm(0.75)*sds), alpha=0.15, color=NA) +
# specify custom colors
scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
labs(x="N", y="", title = "") +
# apply the theme
theme_minimal() +
MyTheme +
# create a facet
ggh4x::facet_nested(factor(netsize, levels = c("5p", "10p"), labels = c("p = 5", "p = 10")) ~ factor(latentvar, levels = c("without LC", "with LC")) + factor(densities, levels=c("sparse", "dense")), scales = "free_y", switch="y") +
scale_x_continuous(breaks=c(50, 2500, 5000, 7500, 10000)) +
ggtitle("SHD") 
See code here.
## Precision figure
pre_rec %>%
filter(grepl("average_precision", metric)) %>%
tidyr::pivot_wider(names_from = metric, values_from=value) %>%
ggplot(aes(x= as.numeric(N), y=average_precision_mean, group = algorithm, colour = algorithm, fill=algorithm)) +
# add line plots
geom_line(aes(group = algorithm)) +
# add scattered points
geom_point(size=1) +
# add interquartile range (IQR)
geom_ribbon(aes(ymin=average_precision_mean+qnorm(0.25)*average_precision_sd, ymax=average_precision_mean+qnorm(0.75)*average_precision_sd), alpha=0.15, color=NA) +
# specify custom colors
scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
# apply the theme
theme_minimal() +
MyTheme +
# create a facet
ggh4x::facet_nested(factor(netsize, levels = c("5p", "10p"), labels=c("p = 5", "p = 10")) ~ factor(latentvar, levels = c("without LC", "with LC")) + factor(densities, levels=c("sparse", "dense")), switch="y") +
scale_x_continuous(breaks=c(50, 2500, 5000, 7500, 10000)) +
labs(title = "Precision", x = "", y = "")
See code here.
## Recall figure
pre_rec %>%
filter(grepl("average_recall", metric)) %>%
tidyr::pivot_wider(names_from = metric, values_from=value) %>%
ggplot(aes(x= as.numeric(N), y=average_recall_mean, group = algorithm, colour = algorithm, fill= algorithm)) +
# add line plots
geom_line(aes(group = algorithm)) +
# add scattered points
geom_point(size=1) +
# add interquartile range (IQR)
geom_ribbon(aes(ymin=average_recall_mean+qnorm(0.25)*average_recall_sd, ymax=average_recall_mean+qnorm(0.75)*average_recall_sd), alpha=0.15, color=NA) +
# specify custom colors
scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
# apply the theme
theme_minimal() +
MyTheme +
# create a facet
ggh4x::facet_nested(factor(netsize, levels = c("5p", "10p"), labels=c("p = 5", "p = 10")) ~ factor(latentvar, levels = c("without LC", "with LC")) + factor(densities, levels=c("sparse", "dense")), switch="y") +
scale_x_continuous(breaks=c(50, 2500, 5000, 7500, 10000)) +
labs(title = "Recall", x = "", y = "")
See code here.
## Uncertainty figure
uncertainties %>%
tidyr::pivot_wider(names_from = statistics, values_from=value) %>%
ggplot(aes(x= as.numeric(N), y=means, group = algorithm, colour = algorithm, fill = algorithm)) +
# add line plots
geom_line(aes(group = algorithm)) +
# add scattered points
geom_point(size=1) +
# add interquartile range (IQR)
geom_ribbon(aes(ymin=means+qnorm(0.25)*sds, ymax=means+qnorm(0.75)*sds), alpha=0.15, color=NA) +
# specify custom colors
scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
scale_fill_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
labs(x="N", y="", title = "") +
# apply the theme
theme_minimal() +
MyTheme +
# create a facet
ggh4x::facet_nested(factor(netsize, levels = c("5p", "10p"), labels=c("p = 5", "p = 10")) ~ factor(latentvar, levels = c("without LC", "with LC")) + factor(densities, levels=c("sparse", "dense")), scales = "free_y", switch="y") +
scale_x_continuous(breaks=c(50, 2500, 5000, 7500, 10000)) +
ggtitle("Uncertainty")
5 Extra: Algorithm Running Time
The running time of the algorithms, presented in log(ms), indicates that FCI generally takes the shortest time, while CCD takes the longest across all conditions.
See code here.
# specify sig. level
alpha <- 0.05
times <- microbenchmark(
ccd_5psparse = ccdKP(df=simdata_5psparse[[10]][[1]], dataType = "continuous", alpha = alpha),
fci_5psparse = fci(list(C = cor(simdata_5psparse[[10]][[1]]), n = 1e1), indepTest=gaussCItest, alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_5psparse[[10]][[1]])),
cci_5psparse = cci(list(C = cor(simdata_5psparse[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_5psparse[[10]][[1]])),
ccd_5pdense = ccdKP(df=simdata_5pdense[[10]][[1]], dataType = "continuous", alpha = alpha),
fci_5pdense = fci(list(C = cor(simdata_5pdense[[10]][[1]]), n = 1e1),indepTest=gaussCItest, alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_5pdense[[10]][[1]])),
cci_5pdense = cci(list(C = cor(simdata_5pdense[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_5pdense[[10]][[1]])),
ccd_10psparse = ccdKP(df=simdata_10psparse[[10]][[1]], dataType = "continuous", alpha = alpha),
fci_10psparse = fci(list(C = cor(simdata_10psparse[[10]][[1]]), n = 1e1),indepTest=gaussCItest, alpha = alpha, doPdsep = TRUE, selectionBias= FALSE, labels = colnames(simdata_10psparse[[10]][[1]])),
cci_10psparse = cci(list(C = cor(simdata_10psparse[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_10psparse[[10]][[1]])),
ccd_10pdense = ccdKP(df=simdata_10pdense[[10]][[1]], dataType = "continuous", alpha = alpha),
fci_10pdense = fci(list(C = cor(simdata_10pdense[[10]][[1]]), n = 1e1),indepTest=gaussCItest, alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_10pdense[[10]][[1]])),
cci_10pdense = cci(list(C = cor(simdata_10pdense[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_10pdense[[10]][[1]])),
ccd_5pLVsparse = ccdKP(df=simdata_5pLVsparse[[10]][[1]], dataType = "continuous", alpha = alpha),
fci_5pLVsparse = fci(list(C = cor(simdata_5pLVsparse[[10]][[1]]), n = 1e1),indepTest=gaussCItest,
alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_5pLVsparse[[10]][[1]])),
cci_5pLVsparse = cci(list(C = cor(simdata_5pLVsparse[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_5pLVsparse[[10]][[1]])),
ccd_5pLVdense = ccdKP(df=simdata_5pLVdense[[10]][[1]], dataType = "continuous", alpha = alpha),
fci_5pLVdense = fci(list(C = cor(simdata_5pLVdense[[10]][[1]]), n = 1e1),indepTest=gaussCItest,
alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_5pLVdense[[10]][[1]])),
cci_5pLVdense = cci(list(C = cor(simdata_5pLVdense[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_5pLVdense[[10]][[1]])),
ccd_10pLVsparse = ccdKP(df=simdata_10pLVsparse[[10]][[1]], dataType = "continuous", alpha = alpha),
fci_10pLVsparse = fci(list(C = cor(simdata_10pLVsparse[[10]][[1]]), n = 1e1),indepTest=gaussCItest,
alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_10pLVsparse[[10]][[1]])),
cci_10pLVsparse = cci(list(C = cor(simdata_10pLVsparse[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_10pLVsparse[[10]][[1]])),
ccd_10pLVdense = ccdKP(df=simdata_10pLVdense[[10]][[1]], dataType = "continuous", alpha = alpha),
fci_10pLVdense = fci(list(C = cor(simdata_10pLVdense[[10]][[1]]), n = 1e1),indepTest=gaussCItest,
alpha = alpha, selectionBias= FALSE, labels = colnames(simdata_10pLVdense[[10]][[1]])),
cci_10pLVdense = cci(list(C = cor(simdata_10pLVdense[[10]][[1]]), n = 1e1), gaussCItest, alpha=alpha, p=ncol(simdata_10pLVdense[[10]][[1]]))
)
times <- times %>%
mutate(algorithm = substr(expr, 1, 3),
condition = stringr::str_split(expr, "_", simplify=T)[,2])
## plot the results
times %>%
ggplot(aes(x=factor(condition, levels= c("5psparse", "5pdense", "10psparse", "10pdense", "5pLVsparse","5pLVdense", "10pLVsparse","10pLVdense")), y = log(time), col= factor(algorithm))) +
geom_boxplot(position = "dodge", outlier.size = 0.8, outlier.alpha = 0.2) + theme_classic() +
# scale_x_discrete(name ="Condition",
# labels=c("", "5p-sparse", "", "","5p-dense","","", "10p-sparse","","","10p-dense","","","5p-LV","","","10p-LV","")) +
scale_colour_manual(values = c("#FF0000", "#00A08A", "#F2AD00"), name= "") +
labs(y = " log(ms)", x = "", title = "Algorithm Running Time", subtitle = "Time in milliseconds (ms)") +
theme(axis.text.x = element_text(face = "bold", family = "Palatino", margin = margin(t = 13), size=10),
legend.position="bottom",
plot.subtitle = element_text(face = "italic", family = "Palatino"),
plot.title = element_text(family = "Palatino", size=15),
axis.title = element_text(family = "Palatino"),
legend.text = element_text(family = "Palatino"))
6 References
Park, Kyuri. 2023. “Discovering cyclic causal models in psychological research.” GitHub repository. https://github.com/KyuriP/Thesis_KP.
Ryan, Oisín, Laura F. Bringmann, and Noémi K. Schuurman. 2022. “The Challenge of Generating Causal Hypotheses Using Network Models.” Structural Equation Modeling: A Multidisciplinary Journal 0 (0): 1–18. https://doi.org/10.1080/10705511.2022.2056039.
7 Session Information
sessionInfo()R version 4.2.3 (2023-03-15)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.0
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] microbenchmark_1.4.9 CCI.KP_0.1.0 MASS_7.3-58.3
[4] ggpubr_0.6.0 ggplot2_3.4.2 qgraph_1.9.4
[7] dplyr_1.1.1 furrr_0.3.1 future_1.32.0
[10] purrr_1.0.1 magrittr_2.0.3 Rgraphviz_2.40.0
[13] graph_1.74.0 BiocGenerics_0.42.0 pcalg_2.7-8
[16] DOT_0.1 rcausal_1.2.1 devtools_2.4.5
[19] usethis_2.1.6 rJava_1.0-6
loaded via a namespace (and not attached):
[1] backports_1.4.1 Hmisc_5.0-1 plyr_1.8.8 igraph_1.4.1
[5] listenv_0.9.0 fastICA_1.2-3 digest_0.6.31 htmltools_0.5.5
[9] fansi_1.0.4 checkmate_2.1.0 memoise_2.0.1 cluster_2.1.4
[13] sfsmisc_1.1-14 remotes_2.4.2 globals_0.16.2 bdsmatrix_1.3-6
[17] prettyunits_1.1.1 ggh4x_0.2.4 jpeg_0.1-10 colorspace_2.1-0
[21] xfun_0.38 callr_3.7.3 crayon_1.5.2 jsonlite_1.8.4
[25] glue_1.6.2 gtable_0.3.3 V8_4.2.2 sjmisc_2.8.9
[29] car_3.1-2 pkgbuild_1.4.0 DEoptimR_1.0-11 ggm_2.5
[33] abind_1.4-5 scales_1.2.1 rstatix_0.7.2 miniUI_0.1.1.1
[37] Rcpp_1.0.10 xtable_1.8-4 htmlTable_2.4.1 clue_0.3-64
[41] foreign_0.8-84 Formula_1.2-5 stats4_4.2.3 profvis_0.3.7
[45] htmlwidgets_1.6.2 lavaan_0.6-15 ellipsis_0.3.2 urlchecker_1.0.1
[49] pkgconfig_2.0.3 farver_2.1.1 nnet_7.3-18 utf8_1.2.3
[53] tidyselect_1.2.0 labeling_0.4.2 rlang_1.1.0 reshape2_1.4.4
[57] later_1.3.0 munsell_0.5.0 tools_4.2.3 cachem_1.0.7
[61] cli_3.6.1 generics_0.1.3 sjlabelled_1.2.0 broom_1.0.4
[65] fdrtool_1.2.17 evaluate_0.20 stringr_1.5.0 fastmap_1.1.1
[69] yaml_2.3.7 processx_3.8.0 knitr_1.42 fs_1.6.1
[73] robustbase_0.95-1 glasso_1.11 pbapply_1.7-0 RBGL_1.72.0
[77] nlme_3.1-162 mime_0.12 compiler_4.2.3 rstudioapi_0.14
[81] curl_5.0.0 png_0.1-8 ggsignif_0.6.4 tibble_3.2.1
[85] pbivnorm_0.6.0 stringi_1.7.12 ps_1.7.4 lattice_0.21-8
[89] Matrix_1.5-4 psych_2.3.3 vctrs_0.6.1 pillar_1.9.0
[93] lifecycle_1.0.3 data.table_1.14.8 cowplot_1.1.1 insight_0.19.1
[97] corpcor_1.6.10 httpuv_1.6.9 R6_2.5.1 promises_1.2.0.1
[101] gridExtra_2.3 parallelly_1.35.0 sessioninfo_1.2.2 codetools_0.2-19
[105] gtools_3.9.4 pkgload_1.3.2 withr_2.5.0 mnormt_2.1.1
[109] parallel_4.2.3 quadprog_1.5-8 rpart_4.1.19 tidyr_1.3.0
[113] rmarkdown_2.21 carData_3.0-5 shiny_1.7.4 base64enc_0.1-3